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Abstract 

This paper presents an algorithm for approximating certain types of 
dynamical systems given by a system of ordinary delay differential equa- 
tions by a Boolean network model. Often Boolean models are much sim- 
pler to understand than complex differential equations models. The mo- 
tivation for this work comes from mathematical systems biology. While 
Boolean mechanisms do not provide information about exact concentra- 
tion rates or time scales, they are often sufficient to capture steady states 
and other key dynamics. Due to their intuitive nature, such models are 
very appealing to researchers in the life sciences. This paper is focused 
on dynamical systems that exhibit bistability and are described by de- 
lay equations. It is shown that if a certain motif including a feedback 
loop is present in the wiring diagram of the system, the Boolean model 
captures the bistability of molecular switches. The method is applied to 
two examples from biology, the lac operon and the phage A lysis/lysogeny 
switch. 



1 Introduction 

Since the discovery of the first gene regulatory network, the lactose metabolism 
network in the bacterium E. Coli by Jacob Monod, [8], such networks have 
been modeled extensively, traditionally with differential equations. But other 
modeling techniques like Boolean networks [f ], stochastic models [15], Petri nets, 
or Bayesian networks [7] have also been used successfully. Work by Kauffman 
in 1969 [9] suggested that gene networks behave like Boolean switching nets 
and therefore Boolean networks are suitable to model them. Using a Boolean 
network model, it has been shown, for instance, that in some cases the network 



topology and the logic of the interactions among the different molecular species 
is sufficient to determine the qualitative dynamics of the network, without taking 
into account the detailed kinetics [f]. In cases where not enough information 
about the kinetic parameters is available to build a detailed continuous model, 
one often still has enough information to build a Boolean network model which 
can provide important information. 

The purpose of this paper is to show for a particular family of continuous 
models that they can be approximated by a Boolean network model that retains 
the key information about model dynamics. Our intent is to demonstrate that 
Boolean models can be used to study a variety of structural and dynamic fea- 
tures of biological and other systems that have traditionally been modeled using 
continuous models. They have the added advantage that they arc intuitive and 
do not require much mathematical background. This makes them particularly 
suitable for use in the life sciences. The focus of this paper is on two features: 
bistability and time delays. We give a general algorithm how to approximate a 
dynamical system given by a system of ordinary delay differential equations by 
a Boolean network model and validate it with two well-known examples from 
systems biology. The relation between discrete and continuous models has been 
examined before, e.g., [6] describes how a logical network can be used to cre- 
ate the analogous differential equations. Other methods allow to systematically 
create a continuous model from a discrete dynamical system [11], or to reduce 
the model space by using a finite state machine [4] . We briefly review the main 
concepts related to Boolean network models. 

1.1 Introduction to Boolean Models 

We use a time discrete deterministic Boolean network with synchronous up- 
date. In a Boolean model, a variable can only be in the "on" (1) or "off" (0) 
state. In systems biology applications, each variable typically corresponds to 
a molecular entity, e.g., the concentration of a gene product such as mRNA 
or proteins. A Boolean network of n variables consists of an update function 
/ = (fi, . . . , /„) : F£ — » F£ , also called a transition function. The dynamics of 
the system are generated deterministically by iteration of the transition func- 
tion /. It can be visualized by its phase space, which shows all possible states 
and their transitions. A cycle in the phase space is called a limit cycle; if it has 
length one, a fixed point. Fixed points of a discrete system are the equivalent 
to steady states of a continuous system. 

In the dependency graph, also referred to as wiring diagram, each variable 
is a vertex and an edge from variable Xi to xj is drawn if Xj shows up in the 
local transition function fj. A directed cycle in the dependency graph is called 
a feedback loop. 
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Example 

Consider the following Boolean network with three variables: 

h = (K)Az 3 ) 

f 2 = (xi A x 3 ) 

/ 3 = ((iiAx 2 )V4 

The phase space of this function is depicted in Fig 1 . The network has a limit 
cycle of length two because the state (1 1) transitions into (Oil), which turns 
back to (1 1). It also includes the fixed point (0 0). Discrete Visualizer of 




Figure 1: Phase space of f = (->xi A 23, x\ A X3, [x\ A £2) V X3) 

Dynamics (DVD tool) [10] was used to create the phase space. It calculates the 
number and length of limit cycles and fixed points, as well as the trajectories 
from any initial state. It also generates a graph of the phase space along with 
the dependency graph. 

2 Approximation by Discrete Models 

We begin by describing several relevant features of biochemical networks that a 
model needs to capture and describe in each case how this is done in the Boolean 
network case. We use lower case letters for continuous variables and upper case 
letters for discrete ones. 

2.0.1 Dilution and Degradation 

It is common in biochemical networks that the concentration of a substance 
X decreases with time because of dilution and degradation. In a differential 
equations model this is usually accounted for with a negative degradation term: 

dx 
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In the discrete model, for simplicity it is assumed that the degradation rate 
for a substance is either vanishingly small and can therefore be ignored, or the 
substance is degraded below the threshold for discretization after two discrete 
timcsteps. 

To model the decrease of concentration of a substance by dilution and degra- 
dation, which is modeled in the continuous case by a negative degradation term 
for x, a new variable X id is introduced. If X \d is "on," then any quantity of 
X that is available has already been reduced once by dilution and degradation, 
and if no new substrate is produced, X will be "off" in the next time step. 

2.0.2 Time Delay 

In gene regulatory networks, time delays are often caused by transcription and 
translation. Dependence of a variable x on the concentration of a substance y 
time r ago can be described with a delay differential equations model 

f t x(t) = ...y(t-T)... 

In a time discrete model, further variables are needed to model this delay. One 
has to choose the length t of a discrete timestep. If the delay r = t, then X is 
chosen to depend on Y\ and Y\ depends on Y. 

fx 

/y = ■■■ 
fvi = Y 

If t > t, additional variables have to be introduced. 

2.0.3 Distinguishing between Low, Medium, and High Concentra- 
tions 

To differentiate between three states for a variable X in a Boolean system, a 
second variable Xhi g h is introduced. By doing so one can make a distinction 
between low, medium, and high concentration of X: if X is on, the concentration 
is a least medium, whereas Xh; g h is only on, if the concentration is high. 

2.0.4 Bistability 

If previous experiments or computations have shown that the system exhibits 
bistable behavior depending on the concentration of a stimulus X, it is nec- 
essary to distinguish in the model between different concentrations of X: the 
concentrations for which there is a unique steady state for the system and for 
which there are multiple possible steady states. Without this distinction, a 
Boolean model could not capture bistability. In the lac operon, a medium al- 
lolactose concentration leads to bistability, therefore the discrete model must 
differentiate between low, medium, and high concentration of allolactose. The 
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Stimulus X Hit* Stimulus X 



\AND /OR 



Figure 2: Necessary motif in dependency graph for bistability 



stimulus X is the parameter that drives a system in and out of the bistable 
region; within the bistable region the previous state determines the resulting 
steady state. The stimulus acts on one or more variables. Usually Xhi g h is part 
of an OR statement in the transition function (if the stimulus is present in high 
concentration, it should overrule everything and turn the system on), whereas 
X is part of an AND statement (under its presence other factors might prevent 
the system from turning on) . The state of Z (this could be a set of variables or 
a single variable) indicates whether the system was in a preinduced state before 
the stimulus concentration was changed or not, depending on the input from Z: 

Y = (Z A X) V X high . 

To assure that Z reflects whether the system was preinduced, Z needs to receive 
input from Y. This leads to a system that can be reduced to a system as the 
one shown in Fig. 2. As long as the dependency graph of a Boolean model can 
be reduced to Fig. 2, it is guaranteed, that the discrete network captures the 
bistability of the system correctly 



2.1 Example of Building Boolean Model 

We demonstrate the above method on the following generic delay differential 
equations model: 

V(t-T)-1 (1) 

x(t - r) - 1. 

The steady state of this system is (x,y) = (1,1). When translating 1 into a 
discrete model, the first step is to generate Boolean equations for x and y, 

fx = y 

fy = X. 

To account for the degradation term —1 in the equation for y, a new variable 
j/oid is introduced as outlined in section 2.0.1, and, similarly, a: id. We use x\ 



dx 
~dt 
dy 
dt 
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and yi for the delay r , as described in section 2.0.2. This leads to the following 
system 

fx = Vi V (x A ^x old ) 

/x ld = "'2/1 • 

/„ = xi V (y A ->y oM ) (2) 

/hi = y 



Part of the phase space of the network is depicted in Fig. 3. It clearly shows 
(x,y) = (1,1) = (1, — , 1, _, _) as its fixed point, just as expected from the solution 
of the continuous model. In the continuous system values for x and y could be 




Figure 3: Phase space of network 2 



negative, this is represented in the fixed point (x,y) = (0,0) = (0, __, 0,_, _) of 
the discrete system. 

In the next two sections we apply the above methods to two well-studied 
examples of gene regulatory networks, for which there exist delay differential 
equations models. 

3 The Lac Operon 

The lac operon is required for the transport and metabolism of lactose in Es- 
cherichia coli. It has been studied extensively, and Novick [12] showed in the 
1950s that bistability is observed with artificial inducers. He observed that 



6 



"prcinduced" bacteria, that is bacteria grown in a high concentration of inducer, 
are able to maintain a high internal inducer concentration if subsequently grown 
in a low external inducer concentration. Novick, [12] calls this the "preinduction 
effect" . He observed, however, that if preinduced bacteria are transferred to a 
medium with no inducer, enzyme synthesis ceases immediately. The minimum 
concentration in which the high rate of synthesis of a preinduced culture can be 
maintained is called the maintenance concentration. 

Novick's experiments have led to various mathematical models of the lac 
operon whose steady state solutions show bistability, for example Yildirim [17], 
Wong[16], Boer [5]. 

We first give a brief overview of the functionality of the lac operon depicted 
in Fig. 4. In the absence of glucose for cellular metabolism, extracellular lactose 
is transported into the cell, cither actively by permease or passively through dif- 
fusion. Inside the cell, /3-galactosidase breaks up lactose into glucose, galactose, 
and allolactose. Allolactose binds to the repressor, which is usually bound to 
the operator region where it inhibits the transcription process, therefore the 
transcription process can proceed. RNA polymerase initiates transcription of 
the structural genes to produce mRNA, which is then translated into proteins, 
permease, and /3-galactosidase. 



Ibc operon 



DMA ■ | 

n n I lacZ I lac Y 
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3 , polymerase 
mRNA . •;: mRNA 
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Protein 



Allolactose ' Inactive 
(inducer) Qr repressor 



p-Galactosidase 



Transacetylase 



(b) Lactose present, repressor inactive, operon on 



Figure 4: Schematic representation of the lactose operon regulatory system from 
dnainfo.wikispaces.com 



The system contains a positive feedback loop, with an increase in the con- 
centration of allolactose, mRNA concentration increases, and therefore more (3- 
galactosidase and permease are produced, which in turn leads to more internal 
lactose and allolactose. Positive (or negative) feedback loops are very common 
in biological systems because they allow for a rapid increase in concentration. 
In the case of the lac operon, the enzyme induction is an "all-or-none phe- 
nomenon" , with the first permease molecule extracellular lactose is transported 
into the cell, increasing the inducer concentration, and therefore increasing the 
probability of the appearance of another permease molecule, [12]. 
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3.1 Boolean Model of Lac Operon 

We derive a Boolean network for the gene regulatory network from the con- 
tinuous model developed in [17] consisting of 5 equations. All its parameters 
were estimated from the biological literature and the equations were numerically 
solved for their steady states. The delay differential equations model predicts 
a bistable region for a medium external inducer concentration, which is in ac- 
cordance to what has been observed experimentally by [12] and [3]. As in the 
continuous model, we take the following 5 variables for the Boolean model into 
account: messenger RNA M, /3-galactosidase B, allolactose A, permease P, and 
internal lactose L. 

We will show how to derive the Boolean equation for messenger RNA, the 
process for the other variables is similar and can be found in the appendix. In 
[17], the rate of change of mRNA is given by the following equation: 

d 1 4- Ki (p-^TM A )n 

— Mit) =a,M^r, tt~, ; — % V 1 o — 7m -M 

dt w m K+Ki(e-^A TM ) n ' 

=aM K + Kl {e-^ A {t-T M )r +T0 ~ (7M + M)M ' 

where n is the number of molecules of allolactose required to inactivate the 
repressor. 4zM{t) depends on A at time (t — tm) and —(7m + fi)M models the 
loss caused by dilution and degradation. 

Under the presence of allolactose A, j^M is non-negative, so the Boolean 
equation for mRNA is 

Jm = XA T , 

where A T describes the allolactose concentration A time tm ago. 

If mRNA is present and no new mRNA is produced in the next time step, the 
concentration will decrease according to the degradation rate. To capture this 
in the Boolean model, we introduce the artificial variable M id as described in 
section 2.0.1. M Q id is "on", if allolactose is absent, because then no new mRNA 
was produced. If no mRNA was produced for 2 time steps, its concentration is 
too low and we consider M to be "off" . This results in the following equations: 

fu = xa t V {x M A -^Moia) 
/m oW = ^X At . 

mRNA depends on A time r ago, so a new variable A\ is introduced, as 
described in 2.0.2 and we set 

hi = XAl V (X M A -^Moic) 
/m oW = ^XAl 
f A = ..- 

Iai = x A - 
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We discretize the external inducer concentration L e to be "on" , if it is above 
the minimal maintenance concentration. The high external inducer concentra- 
tion i ehigh is "on", if the concentration is at least (*), as shown in Fig 5. 



1 

i 



Concentration 

External Inducer 

Figure 5: Discretization of External Inducer 

Using the method described above, accounting for degradation and delays, 
we derive the following Boolean model: 

/A old = V ~^X L 

Jl = (xp A x L J V x L . 

high 

V [x L A ->(af£ old A {x P V x B ))] 
fp = x M 3 V (arp A ->x PtilA ) 

fPald = ^ X M3 

/i e = %L e V X L 

e high 

foe, . . = X L ,. . . ■ 

We use the DVD simulation software, see section 1.1, to calculate that the 
system has 8 fixed points and no limit cycles. 

As expected, a low external inducer concentration (last two variables are 
set to 0) drives the system to fixed points corresponding to the operon in the 
"off" state, points 1 — 3 in table 1. If the concentration of the artificial external 
inducer is medium, but not high, the system results in the fixed points 4 — 7 
in table 1. With high external inducer concentration, the system settles in the 
remaining fixed point which corresponds to an induced operon. 

For a medium inducer concentration, fixed point 4 and 5 are biologically not 
meaningful, because they correspond to states in which only the internal lactose 




Im = x A i V (x M A ->x McM ) 
/a/i = xm 
/a/2 = xmi 

/a/3 = X M2 
/A/oid = ^XAl 
/fl = X M2 V X B 

!a = {xp A xl) V (xl A i£ e ) 

V (x A A ^a;A old ) 
/ai = x A 
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Table 1 : Fixed Points of Lac Operon 



or /3-galactosidase concentration is present but no other substance. Fixed point 
6 corresponds to the "off" operon, fixed point 7 to the "on" operon. 

To show bistability, we analyze the behavior of the fixed points as we change 
the concentration rates of the external inducer. If we start in a state corre- 
sponding to fixed point 1 and increase the inducer concentration to a medium 
concentration, the system settles down in fixed point 6, corresponding to the 
"off" operon. If we start with state 8, the system settles down in fixed point 
7, corresponding to the "on" operon. This is exactly what we expect from the 
solution of the delay differential equation: under a medium concentration of 
external inducer, the steady state depends on whether the cell was preinduced 
or not. 

4 Lambda Phage 

The virus lambda phage (phage A) is a bacteriophage that infects Escherichia 
coli After injecting its DNA into the host, the phage can enter the lytic pathway 
where it alters the host DNA to produce phage particles and then lyses the host 
cell, or it can enter the lysogenic pathway, where it is duplicated with every 
cell division and harmless until the cell is under stress, then it enters its lytic 
pathway. A schematic representation of the phage A switch is shown in Fig. 6. 
Interestingly, the lysogenic state is more stable than the genome itself [2]. A 
comprehensive explanation of the lambda phage switch can be found in [13]. 

In [14] a delay differential equations model for the switch between lysogenic 
and lytic state of phage lambda is presented. This continuous model consists 
of four equations for the rate of change of the concentrations of cl and cro 
messenger RNA molecules, cl and Cro monomers and of two equations for the 
dimer concentrations of [CI2] and [CV02]. 

Bistable behavior depends on the degradation rate of cl, j c i. Numerical 
calculations show that the range for j c j in which the systems has two fixed 
points, one corresponding to lysogeny, the other to lysis, is O.Omin -1 < 7 c j < 
0.35min _1 . If the degradation rate of cl is zero, the phage enters the lysogenic 
pathway, if it is above 0.35min _1 , the lytic pathway. For a medium degradation 
rate of cl, phage lambda can enter the lytic or lysogenic pathway, so one needs 
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to distinguish between three different degradation rates of 7d. Note that for 
the lambda phage a degradation rate drives the system in the bistable region, 
whereas for the lac operon an external inducer has that role. To distinguish 
between three states of 7 c j in the Boolean model, the variable 7 c i M is intro- 
duced. To avoid wrong input, e.g., 7 c i = and 7 c i M h = 1, 7ci is turned on, 
whenever j c i hzgh = 1. 

The delay caused by transcription is only 0.005 min -1 , whereas the delays 
caused by the translation of the monomers cl and Cro are 0.06 and 0.24 min -1 . 
Therefore delay caused by transcription is neglected in the model. 

For a high degradation rate, no cl monomers are left after 1 timestep. We 
assume, that a small concentration of monomers is generated also if the Boolean 
expression that generates cl is zero. As a consequence, if the degradation rate 
for cl is zero, we assume that sufficient cl monomers are produced to set Xci = 1 
after one time step. 
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This results in the following Boolean model. 

/mci = x ci T v x Cro T V (x McI A a;M c i old ) 

/M c i old = ^C/t A ICro T 

fu CIO = x ci^ V (zm cio A x Mcloold ) 

/Afcroold = ^C/t 

x M crol = x M cro 

X Mo,o 2 = X Mara 1 

fci T = TrfT V (7ci Wgh A (a; MoIl V (x C / T A (x C / Told V Tci)))) 

/C/Toid = X ~M~[ 1 V 7cl high 
/Oo T = x M cr02 V (l(7ro T A x Cro Told ) 
fcro ToM — X M 0I02 

Ad = 7cl V 7cl high 
/Tell — 7cl 
/ 7d hish — 7clhi g h ■ 

DVD simulation software [10] is used to calculate the fixed points and gen- 
erate the dynamics in the bistable region. If the degradation rate is O.Omin -1 , 
in the model the last three variables are off, the system results in the fixed point 
[10101001001000] which corresponds to [M c jCI], a high cl and low Cro concen- 
tration which means that the system is in the lysogenic state. If the degradation 
rate is medium, in the model the third last variable is on and the last two vari- 
ables are off, there are two fixed points, [01010110110110] and [10101001001110], 
corresponding to [M cro Cro] and [M c jCI], respectively, lysogenic and lytic path- 
way. With a high degradation rate, in the model only the last variable of the 
last three variables is on, the system settles in the fixed point [01010110110111], 
representing [M cro Cro]. This means that the system is in the lytic state. All 
four fixed points are in accordance with the results found numerically and the 
switching behavior observed experimentally. 

To investigate the bistable region, we start the model with an initialization 
that represents the lysogenic state. Increasing the degradation rate to a medium 
level results in the fixed point [10101001001110], which is still the lysogenic: state. 
Starting the model from a lytic state and decreasing the rate results in the fixed 
point [01010110110110], representing the lytic state. Decreasing the rate even 
further finally results in the fixed point corresponding to the lysogenic state. 

5 Discussion 

The results presented in this paper show that biochemical networks that exhibit 
bistability can be modeled successfully using a Boolean network model, incorpo- 
rating delays for variables, as needed. This was done by showing that continuous 



12 



delay-differential equations models can be approximated by Boolean networks. 
The method presented here is quite general and could be applied to other types 
of biological networks. The examples show that simple Boolean models are able 
to capture steady states and complicated dynamics like hysteresis. 

Boolean network models have the drawback that they do not give rise to 
exact concentration rates of the steady states because their discretization is too 
coarse grained. Discrete models with more than just two states, so called multi 
state models, might be suitable to give enough quantitative information about 
concentration rates while they are intuitive enough and easy to use for a wide 
range of scientists. 
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A Lac Operon 

A.l Dilution and Degradation 

Since the degradation rates 7m, 7p, and ja are close to 0.5mm" 1 , we assume 
that mRNA, permease, and allolactose are degraded after 2 times steps. 71, and 
7s are very small and will be neglected in our model. 

A. 2 /3-galactosidase 

For the /3-galactosidase enzyme the equation in the continuous model is 

^ = a B e-» TB M TB -j B B 

= a B e-^ TB M TB -( lB + n)B. 

Messenger RNA is translated into /3-galactosidase which takes time t b , so f B 
depends on M T , the mRNA concentration time t b ago. 

Since the degradation rate j B is low, the Boolean model neglects the decrease 
due to dilution and degradation and we model /3-galactosidase with the single 
equation 



Jb = x Mt V x B . 



A. 3 Allolactose A 



dA D L A 
dt K L + L Ka + A 
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Allolactosc is gained by conversion of lactose and reduced by the loss via con- 
version to glucose and galactose, both mediated by /3-galactosidase. Like for 
mRNA, we use an extra variable A id to capture the loss of allolactose due to 
dilution and degradation. 

Ia = (xb Ax L )\/ (x A A^x Aald ) (1) 

/Aojd = ^ X B A X L ) 

= ~^x B V ~^x L 

Notice that 1 does not depend on any delayed variables. 
A. 4 Internal Lactose L 

For internal lactose, Yildirim's model suggests the following equation 
dL L e L 

The degradation term 7^ is low enough to be ignored in the discrete model. Lac- 
tose is broken down into glucose, galactose and allolactose by /3-galactosidase, 
and permease actively transports lactose in and out of the cell. Again, since 
permease and /3-galactosidase reduce the internal lactose concentration, we in- 
troduce the extra variable L id to turn L off, if it has not been produced and if 
permease or /3-galactosidase are present to reduce it. 

h = (xp A x Le ) V [x L A ^(x Lold A Op V x B ))] 
/ioii = -^(x P Ax Lc ) 

= ~^Xp V ^Ii c 



A. 5 Permease P 

dP 
dt 



a P e-^ +T ^M Tp+TB -%P 



Messenger RNA is translated into permease, so the permease concentration P 
is directly proportional to the mRNA concentration M at time (rp + tb) ago 
and dilution and degradation reduce permease concentration, which is why P id 
is used in the Boolean model. 

fp = x Mt V (x P A -^x Pold ) 
/Fold = ~> x Mt ■ 
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B Lambda Phage 



Remunerating the equations results in 

/i = % V xio V (xi A Xi) 
f 2 = xs A xio 
h = xi 

h = ^8 V {x 4 A x£) 
h = xb, 
h = x 4 
h = x 6 

h = ^13 V (xu A (x 3 V (x 8 A (Sg V xi^)))) 

/g = V SC14 

/io = x-i V (xio A sll) 

/n = x7 

/12 = x 12 V x 14 

/l3 = £12 
/l4 = X14. 

The Boolean model has the dependency graph shown in Fig. 7. 




Figure 7: Dependency graph of the Boolean model for Lambda Phage 
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